/**********************************************************************
 * 
 *  File: 
 * 
 *  Created:
 *
 *  $Rev: 37 $ 
 *  $Id: Mat4f.hpp 37 2007-08-02 10:21:51Z kchriste $
 *
 *  This file is part of the Sima library.
 *  Copyright (C) 2007 by SINTEF.  All rights reserved.
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  ("GPL") version 2 as published by the Free Software Foundation.
 *  See the file LICENSE.GPL at the root directory of this source
 *  distribution for additional information about the GNU GPL.
 * 
 *  For using Sima with software that can not be combined with the
 *  GNU GPL, please contact SINTEF for aquiring a commercial license
 *  and support.
 *
 *  SINTEF, Pb 124 Blindern, N-0314 Oslo, Norway
 *  http://www.sintef.no
 *********************************************************************/

#ifndef _SIUT_SIMD_STANDARD_MAT4F_HPP_
#define _SIUT_SIMD_STANDARD_MAT4F_HPP_

#include <iostream>
#include <math.h>
#include <cmath>
#include "standard/Vec4f.hpp"
#include "standard/Vec3f.hpp"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif




namespace siut {
  namespace simd {
    /** \brief matrix and vector calculation tool using standard C++ language */

    namespace standard {

      /** \brief Standard: class for doing matrix 4x4 calculations */

      struct Mat4f
      {
	float v_[16];

	Mat4f()
	{

	};

#ifndef NVCC
	Mat4f(const Vec4f &x, const Vec4f &y, const Vec4f &z, const Vec4f &w)
	{
	  v_[0 ] = x.x();
	  v_[1 ] = x.y();
	  v_[2 ] = x.z();
	  v_[3 ] = x.w();
    
	  v_[4 ] = y.x();
	  v_[5 ] = y.y();
	  v_[6 ] = y.z();
	  v_[7 ] = y.w();

	  v_[8 ] = z.x();
	  v_[9 ] = z.y();
	  v_[10] = z.z();
	  v_[11] = z.w();
    
	  v_[12] = w.x();
	  v_[13] = w.y();
	  v_[14] = w.z();
	  v_[15] = w.w();
	}
  
	//by columns
	Mat4f(float m_11, float m_21, float m_31, float m_41, 
		     float m_12, float m_22, float m_32, float m_42,
		     float m_13, float m_23, float m_33, float m_43,
		     float m_14, float m_24, float m_34, float m_44)
	{
	  v_[0 ] = m_11;
	  v_[1 ] = m_21;
	  v_[2 ] = m_31;
	  v_[3 ] = m_41;
    
	  v_[4 ] = m_12;
	  v_[5 ] = m_22;
	  v_[6 ] = m_32;
	  v_[7 ] = m_42;

	  v_[8 ] = m_13;
	  v_[9 ] = m_23;
	  v_[10] = m_33;
	  v_[11] = m_43;
    
	  v_[12] = m_14;
	  v_[13] = m_24;
	  v_[14] = m_34;
	  v_[15] = m_44;  
	}
  
	Mat4f(const float *pf)
	{
	  v_[0 ] = pf[0 ];
	  v_[1 ] = pf[1 ];
	  v_[2 ] = pf[2 ];
	  v_[3 ] = pf[3 ];
    
	  v_[4 ] = pf[4 ];
	  v_[5 ] = pf[5 ];
	  v_[6 ] = pf[6 ];
	  v_[7 ] = pf[7 ];

	  v_[8 ] = pf[8 ];
	  v_[9 ] = pf[9 ];
	  v_[10] = pf[10];
	  v_[11] = pf[11];
    
	  v_[12] = pf[12];
	  v_[13] = pf[13];
	  v_[14] = pf[14];
	  v_[15] = pf[15];
	}
  
	Mat4f(float tab1[4] , float tab2[4] , float tab3[4] , float tab4[4])
	{
	  v_[0 ] = tab1[0 ];
	  v_[1 ] = tab1[1 ];
	  v_[2 ] = tab1[2 ];
	  v_[3 ] = tab1[3 ];
    
	  v_[4 ] = tab2[0 ];
	  v_[5 ] = tab2[1 ];
	  v_[6 ] = tab2[2 ];
	  v_[7 ] = tab2[3 ];

	  v_[8 ] = tab3[0 ];
	  v_[9 ] = tab3[1 ];
	  v_[10] = tab3[2 ];
	  v_[11] = tab3[3 ];
    
	  v_[12] = tab4[0 ];
	  v_[13] = tab4[1 ];
	  v_[14] = tab4[2 ];
	  v_[15] = tab4[3 ];
	}
  
	Mat4f(float tab[16])
	{
	  v_[0 ] = tab[0 ];
	  v_[1 ] = tab[1 ];
	  v_[2 ] = tab[2 ];
	  v_[3 ] = tab[3 ];
    
	  v_[4 ] = tab[4 ];
	  v_[5 ] = tab[5 ];
	  v_[6 ] = tab[6 ];
	  v_[7 ] = tab[7 ];

	  v_[8 ] = tab[8 ];
	  v_[9 ] = tab[9 ];
	  v_[10] = tab[10];
	  v_[11] = tab[11];
    
	  v_[12] = tab[12];
	  v_[13] = tab[13];
	  v_[14] = tab[14];
	  v_[15] = tab[15];
	}
  
	Mat4f(const Mat4f &v)
	{
	  v_[0 ] = v.v_[0 ];
	  v_[1 ] = v.v_[1 ];
	  v_[2 ] = v.v_[2 ];
	  v_[3 ] = v.v_[3 ];
    
	  v_[4 ] = v.v_[4 ];
	  v_[5 ] = v.v_[5 ];
	  v_[6 ] = v.v_[6 ];
	  v_[7 ] = v.v_[7 ];

	  v_[8 ] = v.v_[8 ];
	  v_[9 ] = v.v_[9 ];
	  v_[10] = v.v_[10];
	  v_[11] = v.v_[11];
    
	  v_[12] = v.v_[12];
	  v_[13] = v.v_[13];
	  v_[14] = v.v_[14];
	  v_[15] = v.v_[15];
	}
  
	Mat4f(const float s)  
	{
	  v_[0 ] = s;
	  v_[1 ] = s;
	  v_[2 ] = s;
	  v_[3 ] = s;
    
	  v_[4 ] = s;
	  v_[5 ] = s;
	  v_[6 ] = s;
	  v_[7 ] = s;
    
	  v_[8 ] = s;
	  v_[9 ] = s;
	  v_[10] = s;
	  v_[11] = s;
    
	  v_[12] = s;
	  v_[13] = s;
	  v_[14] = s;
	  v_[15] = s;
	}   
#endif
	const float* c_ptr() const { return &v_[0]; }

	/** Returns the ite element by columns. (obs. first element har index 0) */
	inline float & operator[](int i) { return v_[i]; }
	/** Returns the ite element by columns. (obs. first element har index 0) */
	inline const float & operator[](int i) const { return v_[i]; }
	
	/** Returns element (i,j), where i is column, and j is row) */
	inline float & operator()(int i, int j) { return v_[(i * 4) + j]; }
	/** Returns element (i,j), where i is column, and j is row) */
	inline const float & operator()(int i, int j) const { return v_[(i * 4) + j]; }

      };


      inline void copyInto(Mat4f &ret, const Mat4f &m)
      {
	ret.v_[0] = m.v_[0];
	ret.v_[1] = m.v_[1];
	ret.v_[2] = m.v_[2];
	ret.v_[3] = m.v_[3];
	ret.v_[4] = m.v_[4];
	ret.v_[5] = m.v_[5];
	ret.v_[6] = m.v_[6];
	ret.v_[7] = m.v_[7];
	ret.v_[8] = m.v_[8];
	ret.v_[9] = m.v_[9];
	ret.v_[10] = m.v_[10];
	ret.v_[11] = m.v_[11];
	ret.v_[12] = m.v_[12];
	ret.v_[13] = m.v_[13];
	ret.v_[14] = m.v_[14];
	ret.v_[15] = m.v_[15];
      }

      inline void copyInto(Mat4f &ret, const float s)
      {
	ret.v_[0] = s;
	ret.v_[1] = s;
	ret.v_[2] = s;
	ret.v_[3] = s;
	ret.v_[4] = s;
	ret.v_[5] = s;
	ret.v_[6] = s;
	ret.v_[7] = s;
	ret.v_[8] = s;
	ret.v_[9] = s;
	ret.v_[10] = s;
	ret.v_[11] = s;
	ret.v_[12] = s;
	ret.v_[13] = s;
	ret.v_[14] = s;
	ret.v_[15] = s;
      }



      /*-----------------------------------------------------------------------*/
      /*    operators +=  -=  *= /=                                            */
      /*-----------------------------------------------------------------------*/
 
      inline void operator += (Mat4f &m, const Mat4f &n)
      {
	m.v_[0 ] += n.v_[0 ];
	m.v_[1 ] += n.v_[1 ];
	m.v_[2 ] += n.v_[2 ];
	m.v_[3 ] += n.v_[3 ];
    
	m.v_[4 ] += n.v_[4 ];
	m.v_[5 ] += n.v_[5 ];
	m.v_[6 ] += n.v_[6 ];
	m.v_[7 ] += n.v_[7 ];
    
	m.v_[8 ] += n.v_[8 ];
	m.v_[9 ] += n.v_[9 ];
	m.v_[10] += n.v_[10];
	m.v_[11] += n.v_[11];
    
	m.v_[12] += n.v_[12];
	m.v_[13] += n.v_[13];
	m.v_[14] += n.v_[14];
	m.v_[15] += n.v_[15];
      }
      inline void operator += (Mat4f &m, const float s)
      {
	m.v_[0 ] += s;
	m.v_[1 ] += s;
	m.v_[2 ] += s;
	m.v_[3 ] += s;
    
	m.v_[4 ] += s;
	m.v_[5 ] += s;
	m.v_[6 ] += s;
	m.v_[7 ] += s;
    
	m.v_[8 ] += s;
	m.v_[9 ] += s;
	m.v_[10] += s;
	m.v_[11] += s;
    
	m.v_[12] += s;
	m.v_[13] += s;
	m.v_[14] += s;
	m.v_[15] += s;
      }

      inline void operator -= (Mat4f &m, const Mat4f &n)
      {
	m.v_[0 ] -= n.v_[0 ];
	m.v_[1 ] -= n.v_[1 ];
	m.v_[2 ] -= n.v_[2 ];
	m.v_[3 ] -= n.v_[3 ];
    
	m.v_[4 ] -= n.v_[4 ];
	m.v_[5 ] -= n.v_[5 ];
	m.v_[6 ] -= n.v_[6 ];
	m.v_[7 ] -= n.v_[7 ];
    
	m.v_[8 ] -= n.v_[8 ];
	m.v_[9 ] -= n.v_[9 ];
	m.v_[10] -= n.v_[10];
	m.v_[11] -= n.v_[11];
    
	m.v_[12] -= n.v_[12];
	m.v_[13] -= n.v_[13];
	m.v_[14] -= n.v_[14];
	m.v_[15] -= n.v_[15];
      } 
      inline void operator -= (Mat4f &m, const float s)
      {
	m.v_[0 ] -= s;
	m.v_[1 ] -= s;
	m.v_[2 ] -= s;
	m.v_[3 ] -= s;
    
	m.v_[4 ] -= s;
	m.v_[5 ] -= s;
	m.v_[6 ] -= s;
	m.v_[7 ] -= s;
    
	m.v_[8 ] -= s;
	m.v_[9 ] -= s;
	m.v_[10] -= s;
	m.v_[11] -= s;
    
	m.v_[12] -= s;
	m.v_[13] -= s;
	m.v_[14] -= s;
	m.v_[15] -= s;
      }
  
      inline void operator *= (Mat4f &m, const Mat4f &n)
      {  
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif

	m.v_[ 0] = res.v_[0]*n.v_[ 0] + res.v_[4]*n.v_[ 1] + res.v_[ 8]*n.v_[ 2] + res.v_[12]*n.v_[ 3];
	m.v_[ 1] = res.v_[1]*n.v_[ 0] + res.v_[5]*n.v_[ 1] + res.v_[ 9]*n.v_[ 2] + res.v_[13]*n.v_[ 3];
	m.v_[ 2] = res.v_[2]*n.v_[ 0] + res.v_[6]*n.v_[ 1] + res.v_[10]*n.v_[ 2] + res.v_[14]*n.v_[ 3];
	m.v_[ 3] = res.v_[3]*n.v_[ 0] + res.v_[7]*n.v_[ 1] + res.v_[11]*n.v_[ 2] + res.v_[15]*n.v_[ 3];

	m.v_[ 4] = res.v_[0]*n.v_[ 4] + res.v_[4]*n.v_[ 5] + res.v_[ 8]*n.v_[ 6] + res.v_[12]*n.v_[ 7];
	m.v_[ 5] = res.v_[1]*n.v_[ 4] + res.v_[5]*n.v_[ 5] + res.v_[ 9]*n.v_[ 6] + res.v_[13]*n.v_[ 7];
	m.v_[ 6] = res.v_[2]*n.v_[ 4] + res.v_[6]*n.v_[ 5] + res.v_[10]*n.v_[ 6] + res.v_[14]*n.v_[ 7]; 
	m.v_[ 7] = res.v_[3]*n.v_[ 4] + res.v_[7]*n.v_[ 5] + res.v_[11]*n.v_[ 6] + res.v_[15]*n.v_[ 7];

	m.v_[ 8] = res.v_[0]*n.v_[ 8] + res.v_[4]*n.v_[ 9] + res.v_[ 8]*n.v_[10] + res.v_[12]*n.v_[11];
	m.v_[ 9] = res.v_[1]*n.v_[ 8] + res.v_[5]*n.v_[ 9] + res.v_[ 9]*n.v_[10] + res.v_[13]*n.v_[11];
	m.v_[10] = res.v_[2]*n.v_[ 8] + res.v_[6]*n.v_[ 9] + res.v_[10]*n.v_[10] + res.v_[14]*n.v_[11]; 
	m.v_[11] = res.v_[3]*n.v_[ 8] + res.v_[7]*n.v_[ 9] + res.v_[11]*n.v_[10] + res.v_[15]*n.v_[11];
    
	m.v_[12] = res.v_[0]*n.v_[12] + res.v_[4]*n.v_[13] + res.v_[ 8]*n.v_[14] + res.v_[12]*n.v_[15];
	m.v_[13] = res.v_[1]*n.v_[12] + res.v_[5]*n.v_[13] + res.v_[ 9]*n.v_[14] + res.v_[13]*n.v_[15];
	m.v_[14] = res.v_[2]*n.v_[12] + res.v_[6]*n.v_[13] + res.v_[10]*n.v_[14] + res.v_[14]*n.v_[15]; 
	m.v_[15] = res.v_[3]*n.v_[12] + res.v_[7]*n.v_[13] + res.v_[11]*n.v_[14] + res.v_[15]*n.v_[15];

      }
      inline void operator *= (Mat4f &m, const float s)
      {
	m.v_[0 ] *= s;
	m.v_[1 ] *= s;
	m.v_[2 ] *= s;
	m.v_[3 ] *= s;
    
	m.v_[4 ] *= s;
	m.v_[5 ] *= s;
	m.v_[6 ] *= s;
	m.v_[7 ] *= s;
    
	m.v_[8 ] *= s;
	m.v_[9 ] *= s;
	m.v_[10] *= s;
	m.v_[11] *= s;
    
	m.v_[12] *= s;
	m.v_[13] *= s;
	m.v_[14] *= s;
	m.v_[15] *= s;
      }
  
      inline void operator /= (Mat4f &m, const float s)
      {
	m.v_[0 ] /= s;
	m.v_[1 ] /= s;
	m.v_[2 ] /= s;
	m.v_[3 ] /= s;
    
	m.v_[4 ] /= s;
	m.v_[5 ] /= s;
	m.v_[6 ] /= s;
	m.v_[7 ] /= s;
    
	m.v_[8 ] /= s;
	m.v_[9 ] /= s;
	m.v_[10] /= s;
	m.v_[11] /= s;
    
	m.v_[12] /= s;
	m.v_[13] /= s;
	m.v_[14] /= s;
	m.v_[15] /= s;
      }


      /*-----------------------------------------------------------------------*/
      /*    operators + - * /  with friend                                     */
      /*-----------------------------------------------------------------------*/
  
      inline Mat4f operator + (float s, const Mat4f &m)
      {   
	Mat4f res(s);
	res += m;
	return res;
      }
  
      inline Mat4f operator + (const Mat4f &m, float s)
      {
	return s + m;
      }

      inline Mat4f operator + (const Mat4f &m, const Mat4f &n)
      {
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif

	res += n;
	return res;
      }

      inline Mat4f operator - (float s, const Mat4f &m)
      {   
#ifndef NVCC
	Mat4f res(s);
#else
	Mat4f res;
	copyInto(res, s);
#endif
	res -= m;
	return res;
      }
      
      inline Mat4f operator - (const Mat4f &m, float s)
      {	
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif
	res -= s;
	return res;
      }  

      inline Mat4f operator - (const Mat4f &m, const Mat4f &n)
      {
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif

	res -= n;
	return res;
      }
      
      inline Mat4f operator * (float s, const Mat4f &m)
      {   
#ifndef NVCC
	Mat4f res(m);	
#else
	Mat4f res;
	copyInto(res, m);
#endif
	res *= s;
	return res;
      }

        
      inline Mat4f operator * (const Mat4f &m, float s)
      {
	return s * m;
      }  
      
      inline Vec4f operator * (const Vec4f &v, const Mat4f &m)
      {
	Vec4f res;

	res[0] = v[0]*m.v_[ 0] + v[1]*m.v_[ 1] + v[2]*m.v_[ 2] + v[3]*m.v_[ 3];  
	res[1] = v[0]*m.v_[ 4] + v[1]*m.v_[ 5] + v[2]*m.v_[ 6] + v[3]*m.v_[ 7]; 
	res[2] = v[0]*m.v_[ 8] + v[1]*m.v_[ 9] + v[2]*m.v_[10] + v[3]*m.v_[11]; 
	res[3] = v[0]*m.v_[12] + v[1]*m.v_[13] + v[2]*m.v_[14] + v[3]*m.v_[15]; 

	return res;
      }

      inline Vec4f operator * (const Mat4f &m, const Vec4f &v)
      {
	Vec4f res;
	
	res[0] = v[0] * m.v_[0] + v[1] * m.v_[4] + v[2] * m.v_[8] + v[3] * m.v_[12];
	res[1] = v[0] * m.v_[1] + v[1] * m.v_[5] + v[2] * m.v_[9] + v[3] * m.v_[13];
	res[2] = v[0] * m.v_[2] + v[1] * m.v_[6] + v[2] * m.v_[10] + v[3] * m.v_[14];
	res[3] = v[0] * m.v_[3] + v[1] * m.v_[7] + v[2] * m.v_[11] + v[3] * m.v_[15];

	
	return res;
      }

      inline Mat4f operator * (const Mat4f &m, const Mat4f &n)
      {
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif

	res *= n;
	return res;
      }

      inline Mat4f operator / (const Mat4f &m, const float s)
      {
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif
	res /= s;
	return res;
      }


      /** get - set methods */

      /** \brief returns element (i,j), where i is column, and j is row) */
      inline float& get(Mat4f &m, int i, int j)
      { 
	return m.v_[(i * 4) + j]; 
      }
      /** \brief returns element (i,j), where i is column, and j is row) */
      inline const float& get(const Mat4f &m, int i, int j)
      { 
	return m.v_[(i * 4) + j]; 
      }


      /** \brief returns i-th element */
      inline float& get(Mat4f &m, int i)
      { 
	return m.v_[i]; 
      }
      /** \brief returns i-th element */
      inline const float& get(const Mat4f &m, int i)
      { 
	return m.v_[i];
      }

      inline void set(Mat4f &m, int i, float f)
      {
	m.v_[i] = f;
      }

      inline void set(Mat4f &m, int i, int j, float f)
      {
	m.v_[(i * 4) + j] = f;
      }
  


      /*-----------------------------------------------------------------------*/
      /*    operators == !=                                                    */
      /*-----------------------------------------------------------------------*/

      //operators ==
      inline bool operator == (const Mat4f &m, const Mat4f &u) 
      {    
	for (int i = 0; i < 16; i++)
	  {
	    if (m.v_[i]!=u.v_[i])
	      {
		return false;
	      }
	  }
	return true;

      } 
  
  
      //operator != 
      inline bool operator != (const Mat4f &m, const Mat4f &u) 
      {    
	return !(m == u);
      } 
  

      /*-----------------------------------------------------------------------*/
      /*    functions abs inverse det rotate translate setM...                 */
      /*-----------------------------------------------------------------------*/
 
      inline void abs(Mat4f &m)
      {
	for(int i = 0; i < 16; i++)
	  {
	    if (m.v_[i] < 0)
	      {
		m.v_[i] = m.v_[i] * (-1);
	      }
	  }
      }
  
      inline float det_4x4 (const Mat4f &m)
      {
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif
	float a = res.v_[10] * res.v_[15] - res.v_[14] * res.v_[11];
	float b = res.v_[6 ] * res.v_[15] - res.v_[14] * res.v_[7 ];
	float c = res.v_[6 ] * res.v_[11] - res.v_[10] * res.v_[7 ];
	float d = res.v_[2 ] * res.v_[15] - res.v_[14] * res.v_[3 ];
	float e = res.v_[2 ] * res.v_[11] - res.v_[10] * res.v_[3 ];
	float f = res.v_[2 ] * res.v_[7 ] - res.v_[6 ] * res.v_[3 ];
    
	return(   a * (res.v_[0 ] * res.v_[5 ] - res.v_[4 ] * res.v_[1 ])
		  + b * (res.v_[8 ] * res.v_[1 ] - res.v_[0 ] * res.v_[9 ])
		  + c * (res.v_[0 ] * res.v_[13] - res.v_[12] * res.v_[1 ])
		  + d * (res.v_[4 ] * res.v_[9 ] - res.v_[8 ] * res.v_[5 ])
		  + e * (res.v_[12] * res.v_[5 ] - res.v_[4 ] * res.v_[13])
		  + f * (res.v_[8 ] * res.v_[13] - res.v_[12] * res.v_[9 ])
		  );

      }


      inline Mat4f inverse(const Mat4f &m)                           //116)
      {
	Mat4f res;
    
	float a = ( get(m, 2,2) * get(m, 3,3) - get(m, 2,3) * get(m, 3,2) );       
	float b = ( get(m, 2,0) * get(m, 3,1) - get(m, 2,1) * get(m, 3,0) );
	float c = ( get(m, 3,2) * get(m, 0,3) - get(m, 0,2) * get(m, 3,3) );
	float d = ( get(m, 3,0) * get(m, 0,1) - get(m, 0,0) * get(m, 3,1) );
	float e = ( get(m, 0,2) * get(m, 1,3) - get(m, 0,3) * get(m, 1,2) );
	float f = ( get(m, 0,0) * get(m, 1,1) - get(m, 0,1) * get(m, 1,0) );
	float g = ( get(m, 1,2) * get(m, 2,3) - get(m, 1,3) * get(m, 2,2) );
	float h = ( get(m, 1,0) * get(m, 2,1) - get(m, 1,1) * get(m, 2,0) );
	float i = ( get(m, 1,2) * get(m, 3,3) - get(m, 1,3) * get(m, 3,2) );          
	float j = ( get(m, 1,0) * get(m, 3,1) - get(m, 1,1) * get(m, 3,0) );           
	float k = ( get(m, 0,2) * get(m, 2,3) - get(m, 0,3) * get(m, 2,2) );           
	float l = ( get(m, 0,0) * get(m, 2,1) - get(m, 0,1) * get(m, 2,0) );          
    
	res.v_[ 0] = m.v_[5] * ( a) + m.v_[ 9] * (-i) + m.v_[13] * ( g);    
	res.v_[ 1] = m.v_[1] * (-a) + m.v_[ 9] * (-c) + m.v_[13] * (-k);    
	res.v_[ 2] = m.v_[1] * ( i) + m.v_[ 5] * ( c) + m.v_[13] * ( e);    
	res.v_[ 3] = m.v_[1] * (-g) + m.v_[ 5] * ( k) + m.v_[ 9] * (-e);
    
	res.v_[ 4] = m.v_[4] * (-a) + m.v_[ 8] * ( i) + m.v_[12] * (-g);    
	res.v_[ 5] = m.v_[0] * ( a) + m.v_[ 8] * ( c) + m.v_[12] * ( k);    
	res.v_[ 6] = m.v_[0] * (-i) + m.v_[ 4] * (-c) + m.v_[12] * (-e);    
	res.v_[ 7] = m.v_[0] * ( g) + m.v_[ 4] * (-k) + m.v_[ 8] * ( e);
    
	res.v_[ 8] = m.v_[7] * ( b) + m.v_[11] * (-j) + m.v_[15] * ( h);    
	res.v_[ 9] = m.v_[3] * (-b) + m.v_[11] * (-d) + m.v_[15] * (-l);    
	res.v_[10] = m.v_[3] * ( j) + m.v_[ 7] * ( d) + m.v_[15] * ( f);    
	res.v_[11] = m.v_[3] * (-h) + m.v_[ 7] * ( l) + m.v_[11] * (-f);     
    
	res.v_[12] = m.v_[6] * (-b) + m.v_[10] * ( j) + m.v_[14] * (-h);    
	res.v_[13] = m.v_[2] * ( b) + m.v_[10] * ( d) + m.v_[14] * ( l);    
	res.v_[14] = m.v_[2] * (-j) + m.v_[ 6] * (-d) + m.v_[14] * (-f);    
	res.v_[15] = m.v_[2] * ( h) + m.v_[ 6] * (-l) + m.v_[10] * ( f);          
    
	return res * (float)((1.f)/(det_4x4(m)));
      }
  
      inline float det_3x3(const Mat4f &m) 
      {
#ifndef NVCC
	Mat4f res(m);
#else
	Mat4f res;
	copyInto(res, m);
#endif
	return res.v_[15] *  (   res.v_[0] * (res.v_[5] * res.v_[10] - res.v_[9] * res.v_[ 6]) + 
				 res.v_[4] * (res.v_[9] * res.v_[ 2] - res.v_[1] * res.v_[10]) + 
				 res.v_[8] * (res.v_[1] * res.v_[ 6] - res.v_[5] * res.v_[ 2])
				 );
      }
  
      inline Mat4f inverse_3x3(const Mat4f &m)                           //116)
      {

	Mat4f res;
    
	float a = ( get(m, 2,2) * get(m, 3,3) );       
	float b = ( get(m, 2,0) * get(m, 3,1) - get(m, 2,1) * get(m, 3,0) );
	float c = ( - get(m, 0,2) * get(m, 3,3) );
	float d = ( get(m, 3,0) * get(m, 0,1) - get(m, 0,0) * get(m, 3,1) );
	float f = ( get(m, 0,0) * get(m, 1,1) - get(m, 0,1) * get(m, 1,0) );
	float h = ( get(m, 1,0) * get(m, 2,1) - get(m, 1,1) * get(m, 2,0) );
	float i = ( get(m, 1,2) * get(m, 3,3) );           
	float j = ( get(m, 1,0) * get(m, 3,1) - get(m, 1,1) * get(m, 3,0) );          
	float l = ( get(m, 0,0) * get(m, 2,1) - get(m, 0,1) * get(m, 2,0) );         
    
	res.v_[ 0] = m.v_[5] * ( a) + m.v_[ 9] * (-i) ;    
	res.v_[ 1] = m.v_[1] * (-a) + m.v_[ 9] * (-c) ;    
	res.v_[ 2] = m.v_[1] * ( i) + m.v_[ 5] * ( c) ;    
	res.v_[ 3] = 0;
    
	res.v_[ 4] = m.v_[4] * (-a) + m.v_[ 8] * ( i) ;    
	res.v_[ 5] = m.v_[0] * ( a) + m.v_[ 8] * ( c) ;    
	res.v_[ 6] = m.v_[0] * (-i) + m.v_[ 4] * (-c) ;    
	res.v_[ 7] = 0;
    
	res.v_[ 8] = m.v_[15] * ( h);    
	res.v_[ 9] = m.v_[15] * (-l);    
	res.v_[10] = m.v_[15] * ( f);    
	res.v_[11] = m.v_[3] * (-h) + m.v_[ 7] * ( l) + m.v_[11] * (-f);     

	res.v_[12] = m.v_[6] * (-b) + m.v_[10] * ( j) + m.v_[14] * (-h);    
	res.v_[13] = m.v_[2] * ( b) + m.v_[10] * ( d) + m.v_[14] * ( l);    
	res.v_[14] = m.v_[2] * (-j) + m.v_[ 6] * (-d) + m.v_[14] * (-f);    
	res.v_[15] = m.v_[2] * ( h) + m.v_[ 6] * (-l) + m.v_[10] * ( f);          
    
	return res * (float(1.0)/(det_3x3(m)));
      }
  
  
   
  
      inline Mat4f inverse_2 (const Mat4f &m)                  //(272)
      {
	Mat4f res;
	
	set(res, 0, 0, get(m, 1,2)*get(m, 2,3)*get(m, 3,1) - get(m, 1,3)*get(m, 2,2)*get(m, 3,1) 
	    + get(m, 1,3)*get(m, 2,1)*get(m, 3,2) - get(m, 1,1)*get(m, 2,3)*get(m, 3,2) 
	    - get(m, 1,2)*get(m, 2,1)*get(m, 3,3) + get(m, 1,1)*get(m, 2,2)*get(m, 3,3));    
	set(res, 0, 1, get(m, 0,3)*get(m, 2,2)*get(m, 3,1) - get(m, 0,2)*get(m, 2,3)*get(m, 3,1) 
	    - get(m, 0,3)*get(m, 2,1)*get(m, 3,2) + get(m, 0,1)*get(m, 2,3)*get(m, 3,2) 
	    + get(m, 0,2)*get(m, 2,1)*get(m, 3,3) - get(m, 0,1)*get(m, 2,2)*get(m, 3,3));    
	set(res, 0, 2, get(m, 0,2)*get(m, 1,3)*get(m, 3,1) - get(m, 0,3)*get(m, 1,2)*get(m, 3,1) 
	    + get(m, 0,3)*get(m, 1,1)*get(m, 3,2) - get(m, 0,1)*get(m, 1,3)*get(m, 3,2) 
	    - get(m, 0,2)*get(m, 1,1)*get(m, 3,3) + get(m, 0,1)*get(m, 1,2)*get(m, 3,3));    
	set(res, 0, 3, get(m, 0,3)*get(m, 1,2)*get(m, 2,1) - get(m, 0,2)*get(m, 1,3)*get(m, 2,1) 
	  - get(m, 0,3)*get(m, 1,1)*get(m, 2,2) + get(m, 0,1)*get(m, 1,3)*get(m, 2,2) 
	    + get(m, 0,2)*get(m, 1,1)*get(m, 2,3) - get(m, 0,1)*get(m, 1,2)*get(m, 2,3));    
	set(res, 1, 0, get(m, 1,3)*get(m, 2,2)*get(m, 3,0) - get(m, 1,2)*get(m, 2,3)*get(m, 3,0) 
	  - get(m, 1,3)*get(m, 2,0)*get(m, 3,2) + get(m, 1,0)*get(m, 2,3)*get(m, 3,2) 
	    + get(m, 1,2)*get(m, 2,0)*get(m, 3,3) - get(m, 1,0)*get(m, 2,2)*get(m, 3,3));    
	set(res, 1, 1, get(m, 0,2)*get(m, 2,3)*get(m, 3,0) - get(m, 0,3)*get(m, 2,2)*get(m, 3,0) 
	  + get(m, 0,3)*get(m, 2,0)*get(m, 3,2) - get(m, 0,0)*get(m, 2,3)*get(m, 3,2) 
	    - get(m, 0,2)*get(m, 2,0)*get(m, 3,3) + get(m, 0,0)*get(m, 2,2)*get(m, 3,3));    
	set(res, 1, 2, get(m, 0,3)*get(m, 1,2)*get(m, 3,0) - get(m, 0,2)*get(m, 1,3)*get(m, 3,0) 
	  - get(m, 0,3)*get(m, 1,0)*get(m, 3,2) + get(m, 0,0)*get(m, 1,3)*get(m, 3,2) 
	    + get(m, 0,2)*get(m, 1,0)*get(m, 3,3) - get(m, 0,0)*get(m, 1,2)*get(m, 3,3));    
	set(res, 1, 3, get(m, 0,2)*get(m, 1,3)*get(m, 2,0) - get(m, 0,3)*get(m, 1,2)*get(m, 2,0) 
	  + get(m, 0,3)*get(m, 1,0)*get(m, 2,2) - get(m, 0,0)*get(m, 1,3)*get(m, 2,2) 
	    - get(m, 0,2)*get(m, 1,0)*get(m, 2,3) + get(m, 0,0)*get(m, 1,2)*get(m, 2,3));    
	set(res, 2, 0, get(m, 1,1)*get(m, 2,3)*get(m, 3,0) - get(m, 1,3)*get(m, 2,1)*get(m, 3,0) 
	  + get(m, 1,3)*get(m, 2,0)*get(m, 3,1) - get(m, 1,0)*get(m, 2,3)*get(m, 3,1) 
	    - get(m, 1,1)*get(m, 2,0)*get(m, 3,3) + get(m, 1,0)*get(m, 2,1)*get(m, 3,3));    
	set(res, 2, 1, get(m, 0,3)*get(m, 2,1)*get(m, 3,0) - get(m, 0,1)*get(m, 2,3)*get(m, 3,0) 
	  - get(m, 0,3)*get(m, 2,0)*get(m, 3,1) + get(m, 0,0)*get(m, 2,3)*get(m, 3,1) 
	    + get(m, 0,1)*get(m, 2,0)*get(m, 3,3) - get(m, 0,0)*get(m, 2,1)*get(m, 3,3));    
	set(res, 2, 2, get(m, 0,1)*get(m, 1,3)*get(m, 3,0) - get(m, 0,3)*get(m, 1,1)*get(m, 3,0) 
	  + get(m, 0,3)*get(m, 1,0)*get(m, 3,1) - get(m, 0,0)*get(m, 1,3)*get(m, 3,1) 
	    - get(m, 0,1)*get(m, 1,0)*get(m, 3,3) + get(m, 0,0)*get(m, 1,1)*get(m, 3,3));    
	set(res, 2, 3, get(m, 0,3)*get(m, 1,1)*get(m, 2,0) - get(m, 0,1)*get(m, 1,3)*get(m, 2,0) 
	  - get(m, 0,3)*get(m, 1,0)*get(m, 2,1) + get(m, 0,0)*get(m, 1,3)*get(m, 2,1) 
	    + get(m, 0,1)*get(m, 1,0)*get(m, 2,3) - get(m, 0,0)*get(m, 1,1)*get(m, 2,3));    
	set(res, 3, 0, get(m, 1,2)*get(m, 2,1)*get(m, 3,0) - get(m, 1,1)*get(m, 2,2)*get(m, 3,0) 
	  - get(m, 1,2)*get(m, 2,0)*get(m, 3,1) + get(m, 1,0)*get(m, 2,2)*get(m, 3,1) 
	    + get(m, 1,1)*get(m, 2,0)*get(m, 3,2) - get(m, 1,0)*get(m, 2,1)*get(m, 3,2));    
	set(res, 3, 1, get(m, 0,1)*get(m, 2,2)*get(m, 3,0) - get(m, 0,2)*get(m, 2,1)*get(m, 3,0) 
	  + get(m, 0,2)*get(m, 2,0)*get(m, 3,1) - get(m, 0,0)*get(m, 2,2)*get(m, 3,1) 
	    - get(m, 0,1)*get(m, 2,0)*get(m, 3,2) + get(m, 0,0)*get(m, 2,1)*get(m, 3,2));    
	set(res, 3, 2, get(m, 0,2)*get(m, 1,1)*get(m, 3,0) - get(m, 0,1)*get(m, 1,2)*get(m, 3,0) 
	  - get(m, 0,2)*get(m, 1,0)*get(m, 3,1) + get(m, 0,0)*get(m, 1,2)*get(m, 3,1) 
	    + get(m, 0,1)*get(m, 1,0)*get(m, 3,2) - get(m, 0,0)*get(m, 1,1)*get(m, 3,2));    
	set(res, 3, 3, get(m, 0,1)*get(m, 1,2)*get(m, 2,0) - get(m, 0,2)*get(m, 1,1)*get(m, 2,0) 
	  + get(m, 0,2)*get(m, 1,0)*get(m, 2,1) - get(m, 0,0)*get(m, 1,2)*get(m, 2,1) 
	    - get(m, 0,1)*get(m, 1,0)*get(m, 2,2) + get(m, 0,0)*get(m, 1,1)*get(m, 2,2));    
    
	return res * (1.f/det_4x4(m));
      }
  

      inline void setIdentity(Mat4f &m) 
      {
	m.v_[0 ] = 1;
	m.v_[1 ] = 0;
	m.v_[2 ] = 0;
	m.v_[3 ] = 0;
    
	m.v_[4 ] = 0;
	m.v_[5 ] = 1;
	m.v_[6 ] = 0;
	m.v_[7 ] = 0;
    
	m.v_[8 ] = 0;
	m.v_[9 ] = 0;
	m.v_[10] = 1;
	m.v_[11] = 0;
    
	m.v_[12] = 0;
	m.v_[13] = 0;
	m.v_[14] = 0;
	m.v_[15] = 1;
      }
  

      //satt inn enere på diagonalen i rotasjonsmatrisene (andré har godkjent)
      inline void rotateX(Mat4f &m, float a)
      {
	
	Mat4f tmp;
	setIdentity(tmp);
	float cosa=cos(a);
	float sina=sin(a);
	tmp.v_[5]=cosa;
	tmp.v_[6]=sina;
	tmp.v_[9]=-sina;
	tmp.v_[10]=cosa;
	m *= tmp;	
      }
  
      inline void rotateY(Mat4f &m, float a)
      {
	Mat4f tmp;
	setIdentity(tmp);
	float cosa=cos(a);
	float sina=sin(a);
	tmp.v_[0]=cosa;
	tmp.v_[2]=-sina;
	tmp.v_[8]=sina;
	tmp.v_[10]=cosa;
	m *= tmp;
      }
  
      inline void rotateZ(Mat4f &m, float a)
      {
	Mat4f tmp;
	setIdentity(tmp);
	float cosa=cos(a);
	float sina=sin(a);
	tmp.v_[0]=cosa;
	tmp.v_[1]=sina;
	tmp.v_[4]=-sina;
	tmp.v_[5]=cosa;
	m *= tmp;
      }
  
      inline void translate(Mat4f &m, float x, float y, float z)
      {
	Mat4f tmp;
	setIdentity(tmp);
	tmp.v_[12]=x;
	tmp.v_[13]=y;
	tmp.v_[14]=z;
	m *= tmp;
      }
  
  
      inline void translate(Mat4f &m, const Vec3f &v)
      {
	translate(m, v[0], v[1], v[2]);
      }
  
     
      inline void setFrustumMatrix(Mat4f &m, float leftV, float rightV, float bottomV, float topV, float nearV, float farV)
      {
	m.v_[0]  = (2.0f * nearV) / (rightV - leftV);
	m.v_[1]  = 0.0;
	m.v_[2]  = 0.0;
	m.v_[3]  = 0.0;
     
	m.v_[4]  = 0.0;
	m.v_[5]  = (2.0f * nearV) / (topV - bottomV);
	m.v_[6]  = 0.0;
	m.v_[7]  = 0.0;
    
	m.v_[8]  = (rightV + leftV) / (rightV - leftV);
	m.v_[9]  = (topV + bottomV) / (topV - bottomV);
	m.v_[10] = -(farV + nearV) / (farV - nearV);
	m.v_[11] = -1.0f;
    
	m.v_[12] = 0.0;
	m.v_[13] = 0.0;
	m.v_[14] = -(2.0f * farV * nearV) / (farV - nearV);
	m.v_[15] = 0.0;
      }
  
      inline void setPerspectiveMatrix(Mat4f &m, float fovy, float aspect, float zNear, float zFar)
      {
	float xMin, xMax, yMin, yMax;
    
	yMax = static_cast<float>(zNear * (tan(fovy * M_PI / 360.0)));
	yMin = -yMax;
	xMin = yMin * aspect;
	xMax = yMax * aspect;
    
	setFrustumMatrix(m, xMin, xMax, yMin, yMax, zNear, zFar);
      }
  
      inline Vec3f rotateVector(Mat4f &m, const Vec3f &v)
      {
	Vec3f vec;
	vec[0] = m.v_[0] * v[0] + m.v_[1] * v[1] + m.v_[ 2] * v[2];
	vec[1] = m.v_[4] * v[0] + m.v_[5] * v[1] + m.v_[ 6] * v[2];
	vec[2] = m.v_[8] * v[0] + m.v_[9] * v[1] + m.v_[10] * v[2];
	return vec;
      }
  
	  // Added by Torbjoern
	  // this is allows one to rotate and translate a Vec3f by a 4x4 matrix
	  inline Vec3f transform(const Mat4f &m, const Vec3f &v)
	  {
		  Vec3f vec;
		  vec[0] = m.v_[0] * v[0] + m.v_[1] * v[1] + m.v_[ 2] * v[2] + m.v_[3];
		  vec[1] = m.v_[4] * v[0] + m.v_[5] * v[1] + m.v_[ 6] * v[2] + m.v_[7];
		  vec[2] = m.v_[8] * v[0] + m.v_[9] * v[1] + m.v_[10] * v[2] + m.v_[11];
		  return vec;
	  }

 
      /*-----------------------------------------------------------------------*/
      /*    functions c1 c2 c3 c4 r1 r2 r3 r4                                  */
      /*-----------------------------------------------------------------------*/
  
    
      inline Vec4f c1(const Mat4f &m)  
      {
	Vec4f tmp(m.v_[0],m.v_[1],m.v_[2],m.v_[3]);
	return tmp;  
      }
  
     
      inline Vec4f c2(const Mat4f &m)  
      {
	Vec4f tmp(m.v_[4],m.v_[5],m.v_[6],m.v_[7]);
	return tmp;  
      }
  
     
      inline Vec4f c3(const Mat4f &m)  
      {
	Vec4f tmp(m.v_[8],m.v_[9],m.v_[10],m.v_[11]);
	return tmp;  
      }
  
      inline Vec4f c4(const Mat4f &m)  
      {        
	Vec4f tmp(m.v_[12],m.v_[13],m.v_[14],m.v_[15]);
	return tmp;  
      }
  
      inline Vec4f r1(const Mat4f &m)  
      {
	Vec4f tmp(m.v_[0],m.v_[4],m.v_[8],m.v_[12]);
	return tmp;  
      }
  
      inline Vec4f r2(const Mat4f &m)  
      {
	Vec4f tmp(m.v_[1],m.v_[5],m.v_[9],m.v_[13]);
	return tmp;  
      }
  
     
      inline Vec4f r3(const Mat4f &m)  
      {
	Vec4f tmp(m.v_[2],m.v_[6],m.v_[10],m.v_[14]);
	return tmp;  
      }
  
     
      inline Vec4f r4(const Mat4f &m)  
      {
	Vec4f tmp(m.v_[3],m.v_[7],m.v_[11],m.v_[15]);
	return tmp;  
      }

      /*-----------------------------------------------------------------------*/
      /*    read  write                                                        */
      /*-----------------------------------------------------------------------*/

      //by colum
      inline void read(Mat4f &m, std::istream& is)
      {
	is >> m.v_[0] >> m.v_[4] >> m.v_[8]  >> m.v_[12]
	   >> m.v_[1] >> m.v_[5] >> m.v_[9]  >> m.v_[13]
	   >> m.v_[2] >> m.v_[6] >> m.v_[10] >> m.v_[14]
	   >> m.v_[3] >> m.v_[7] >> m.v_[11] >> m.v_[15];
      }
  
      inline void write(const Mat4f &m, std::ostream& os) 
      {
	os <<"\n"
	   <<m.v_[0]<<" "<<m.v_[4]<<" "<<m.v_[8] <<" "<<m.v_[12]<<"\n"
	   <<m.v_[1]<<" "<<m.v_[5]<<" "<<m.v_[9] <<" "<<m.v_[13]<<"\n"
	   <<m.v_[2]<<" "<<m.v_[6]<<" "<<m.v_[10]<<" "<<m.v_[14]<<"\n"
	   <<m.v_[3]<<" "<<m.v_[7]<<" "<<m.v_[11]<<" "<<m.v_[15];
      }
  
   
  
      /*-----------------------------------------------------------------------*/
      /*    functions max min tr with friend                                   */
      /*-----------------------------------------------------------------------*/  
  
      //func max/min with friend
      inline Mat4f min (const Mat4f &u , const Mat4f &v) 
      { 
	Mat4f res;  

	for (int i=0; i<16; i++)
	  {
	    if(u.v_[i]<v.v_[i])
	      res.v_[i]=u.v_[i];
	    else
	      res.v_[i]=v.v_[i];
	  }

	return res;
      }
      inline Mat4f max (const Mat4f &u , const Mat4f &v) 
      {   
	Mat4f res;  

	for (int i=0; i<16; i++)
	  {
	    if(u.v_[i]>=v.v_[i])
	      res.v_[i]=u.v_[i];
	    else
	      res.v_[i]=v.v_[i];
	  }

	return res;
      }
  
      inline Mat4f transpose (const Mat4f &m)
      {
	Mat4f res(m.v_[0], m.v_[4], m.v_[8 ], m.v_[12],
		  m.v_[1], m.v_[5], m.v_[9 ], m.v_[13],
		  m.v_[2], m.v_[6], m.v_[10], m.v_[14],
		  m.v_[3], m.v_[7], m.v_[11], m.v_[15]);
	return res;
      }

      inline float trace (const Mat4f &m) 
      {
	float res;
	res=m.v_[0] + m.v_[5] + m.v_[10] + m.v_[15];    

	return res;
      }


      inline std::istream& operator>>(std::istream& is, Mat4f& m)
      { 
	read(m, is); 
	return is; 
      }
  

      inline std::ostream& operator<<(std::ostream& os, const Mat4f& m)
      { 
	write(m, os); 
	return os; 
      }  

      /*-----------------------------------------------------------------------*/
      /*    functions perspective-/frustum-/orthographic-matrix()              */
      /*-----------------------------------------------------------------------*/  
  

      inline void perspectiveMat4f(Mat4f &m, float fovY, float aspect, float zNear, float zFar)
      {
	float f = static_cast<float>(1.0/tan(0.5 * fovY));
	float g = f/aspect;
	float zTemp = zNear - zFar;
	float h = (zFar + zNear) / zTemp;
	float j = (2 * zFar * zNear) / zTemp;

	Mat4f temp(g, 0.f, 0.f, 0.f, 0.f, f, 0.f, 0.f, 0.f, 0.f, h, j, 0.f, 0.f, -1.f, 0.f);
	m = temp;
      }

      inline void orthographicMat4f(Mat4f &m, float left, float right, float bot, float top, float lNear, float lFar)
      {
	float trl = right - left;
	float ttb = top - bot;
	float tfn = lFar - lNear;
	float tx = -(right + left) / trl;
	float ty = -(top + bot) / ttb;
	float tz = -(lFar + lNear) / tfn;
	float sx = 2.f / trl;
	float sy = 2.f / ttb;
	float sz = -2.f / tfn;
    
	Mat4f temp(sx, 0.f, 0.f, tx, 0.f, sy, 0.f, ty, 0.f, 0.f, sz, tz, 0.f, 0.f, 0.f, 1.f);
	m = temp;
      }

      inline void frustumMat4f(Mat4f &m, float left, float right, float bot, float top, float lNear, float lFar)
      {
	float trl = right - left;
	float ttb = top - bot;
	float tfn = lFar - lNear;
    
	float a = (right + left) / trl;
	float b = (top + bot) / ttb;
	float c = -(lFar + lNear) / tfn;
	float d = -2.f * lNear / tfn;
	float e = 2.f * lNear / trl;
	float f = 2.f * lNear / ttb;


	Mat4f temp(e, 0.f, a, 0.f, 0.f, f, b, 0.f, 0.f, 0.f, c, d, 0.f, 0.f, -1.f, 0.f);
	m = temp;
      }


    }//end namespaces
  }
}
#endif
